Source code for hysop.numerics.stencil.stencil

# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


"""Finite differences stencil inteface used in various numerical methods.

* :class:`~hysop.numerics.stencil.stencil.Stencil`
* :class:`~hysop.numerics.stencil.stencil.CenteredStencil`
* :class:`~hysop.numerics.stencil_generator.StencilGenerator`
* :class:`~hysop.numerics.stencil_generator.CenteredStencilGenerator`

"""

import hashlib
import numpy as np
import scipy as sp
import sympy as sm
import itertools as it

from hysop.tools.htypes import check_instance, first_not_None, to_tuple
from hysop.tools.sympy_utils import recurse_expression_tree


[docs] class Stencil: """ Stencil used for finite abitrary dimension differences. """ def __init__( self, coeffs, origin, order, dx=sm.Symbol("dx"), factor=1, error=None, delete_zeros=True, ): """ Stencil used for finite differences schemes (n-dimensional). Parameters ---------- coeffs: array_like Coefficients of the stencil. origin: array_like Origin of the stencil from left-most point. order: int Order of the stencil. Other Parameters ---------------- dx: sympy.Symbol or array_like of sympy.Symbol Spatial spacing symbols. factor: dtype or sympy.Expr Common leading factors for all coefficients. error: sympy.Expr Exact error if known. Attributes ---------- coeffs: np.ndarray Coefficients of the stencil. origin: np.ndarray Origin of the stencil from left-most point. order: int Spatial order of the stencil. dim: int Spatial dimension of the stencil. shape: np.ndarray Shape of the stencil. L: np.ndarray Distance from origin to left-most point included in stencil support. R: np.ndarray Distance from origin to right-most point included in stencil support. dx: sympy.Symbol or array_like of sympy.Symbol Spatial spacing symbols. factor: dtype or sympy.Expr Leading factor of all coefficients. error: sympy.Expr Exact error term if known or else None. Raises ------ ValueError Raised if one component of `origin` is negative. See Also -------- :class:`CenteredStencil`: Centered version of a stencil. :class:`StencilGenerator`: Generate Stencil objects. """ coeffs = np.atleast_1d(coeffs) origin = np.atleast_1d(origin) order = np.atleast_1d(order) dx = np.atleast_1d(dx) if (origin < 0).any(): raise ValueError(f"Origin component < 0!\norigin={origin}") self.dx = dx[0] if dx.size == 1 else dx self.error = error self.origin = origin[0] if origin.size == 1 else origin self.order = order[0] if order.size == 1 else order self.factor = factor if delete_zeros: coeffs = self._delete_zeros(coeffs) self.coeffs = coeffs self._update_attributes()
[docs] def reshape(self, new_shape): """Reshape a stencil by adding zeros.""" new_shape = np.asarray(new_shape) shape = np.asarray(self.shape) assert new_shape.ndim == shape.ndim assert (new_shape >= shape).all() assert ((new_shape - shape) % 2 == 0).all() zeros = (new_shape - shape) // 2 slc = tuple(slice(z, z + s, 1) for (z, s) in zip(zeros, shape)) new_origin = zeros + self.origin new_coeffs = np.zeros(shape=new_shape, dtype=self.coeffs.dtype) new_coeffs[slc] = self.coeffs return self.__class__( coeffs=new_coeffs, origin=new_origin, order=self.order, dx=self.dx, factor=self.factor, error=self.error, delete_zeros=False, )
[docs] def has_factor(self): return self.factor != 1
[docs] def non_zero_coefficients(self): return np.sum(self.coeffs != 0)
[docs] def replace_symbols(self, dic): if isinstance(self.factor, sm.Basic): self.factor = self.factor.xreplace(dic) coeffs = self.coeffs.ravel() for i, coeff in enumerate(coeffs): if isinstance(coeff, sm.Basic): coeffs[i] = coeff.xreplace(dic)
[docs] def apply(self, a, out=None, axis=None, symbols=None, iview=Ellipsis): """ Apply this stencil the input array a. Boundaries are considered periodic. """ symbols = first_not_None(symbols, {}) check_instance(a, np.ndarray) check_instance(symbols, dict) adim = a.ndim sdim = self.dim assert sdim <= adim, "Stencil dimension greater than array dimension." assert ( set(symbols.keys()) == self.variables() ), f"Missing symbols {self.variables()-set(symbols.keys())}." out = first_not_None(out, np.empty_like(a[iview])) axis = first_not_None(to_tuple(axis), tuple(range(adim))[-sdim:]) assert len(axis) == sdim assert out.ndim == a.ndim assert out.shape == a[iview].shape out[...] = 0 for shift, coeff in self.items(include_factor=False): if isinstance(coeff, sm.Basic): coeff = coeff.xreplace(symbols) coeff = float(coeff) view = a.view() for i, si in enumerate(shift): view = np.roll(view, shift=-si, axis=axis[i]) out[...] += coeff * view[iview] factor = self.factor if isinstance(self.factor, sm.Basic): factor = self.factor.xreplace(symbols) factor = float(factor) out[...] *= factor return out
[docs] def __call__(self, *args, **kwds): """Alias for apply.""" return self.apply(*args, **kwds)
def _update_attributes(self): self.dim = self.coeffs.ndim self.shape = self.coeffs.shape self.L = np.asarray([self.origin] if np.isscalar(self.origin) else self.origin) self.R = self.shape - self.origin - 1 def _delete_zeros(self, coeffs): dim = coeffs.ndim shape = list(coeffs.shape) dtype = coeffs.dtype if dtype in [np.float16, np.float32, np.float64]: mask = np.isclose(coeffs, 0, atol=2 * np.finfo(dtype).eps) else: mask = coeffs == 0 keep_mask = np.ones_like(coeffs, dtype=bool) for d in range(dim): laccess = [slice(None)] * dim raccess = [slice(None)] * dim laccess[d] = slice(0, 1) raccess[d] = slice(-1, None) laccess = tuple(laccess) raccess = tuple(raccess) ldel = mask[laccess].all() rdel = mask[raccess].all() if ldel: keep_mask[tuple(laccess)] = False if dim == 1: self.origin -= 1 else: self.origin[d] -= 1 shape[d] -= 1 if rdel: shape[d] -= 1 keep_mask[tuple(raccess)] = False coeffs = coeffs[keep_mask].reshape(shape) return coeffs
[docs] def accuracy(self): """ Accuracy of the stencil with bigO notation. Returns ------- accuracy: sympy.Expr Accuracy in bigO notation. """ expr = 0 done = [] if self.dim != 1: for dx, order in zip(self.dx, self.order): if dx not in done: expr += sm.O(dx ** int(order)) done.append(dx) else: expr += sm.O(self.dx ** int(self.order)) return expr
[docs] def variables(self): """ Symbolic variables used by the stencil. Returns ------- vars: set of sympy.Expr """ _vars = set() def op(expr): if isinstance(expr, sm.Symbol): _vars.add(expr) recurse_expression_tree(op, self.factor) for coeff in self.coeffs.flatten(): recurse_expression_tree(op, coeff) return _vars
[docs] def coo_stencil(self): """ Return a 2d stencil as a sparse coordinate matrix. Returns ------- coo_stencil: scipy.sparse.coo_matrix Sparse Coordinate Matrix of the stencil. Raises ------ RuntimeError Raised if stencil is not 2d. """ if self.dim != 2: raise RuntimeError("Stencil is not 2d !") return sp.sparse.coo_matrix(self.coeffs, shape=self.shape, dtype=self.dtype)
[docs] def items(self, include_factor=True): """ Return an (offset,coefficient) iterator iterating on all **non zero** coefficients. Offset is taken from origin. Returns ------- it : itertools.iterator Zipped offset and coefficient iterator. """ factor = self.factor if include_factor else 1 def mapfun(x): offset = x - self.origin value = factor * self.coeffs[x] return (offset, value) iterator = np.ndindex(self.shape) iterator = map(mapfun, iterator) iterator = filter(lambda x: x[1] != 0, iterator) return iterator
[docs] def refactor(self, factor): """ Factorize the stencil by a new numeric or symbolic value. Parameters ---------- factor : :attr:`dtype` or sympy.Expr New factor to be applied. Returns ------- self : :class:`Stencil` """ old_factor = self.factor cfactor = old_factor / factor coeffs = self.coeffs for idx, coeff in np.ndenumerate(coeffs): coeffs[idx] = sm.simplify(coeff * cfactor) self.factor = factor self.coeffs = coeffs return self
[docs] def hash(self, keep_only=16): """ Hash the stencil with its origin, order and coefficients. The hash algorithm used is sha1. By default only the 16 first chars are kept from the generated hash. Parameters ---------- keep_only : int Number of chars kept from the hashed hexedecimal string. Returns ------- hash : string Hexadecimal hash string of size :attr:`keep_only`. """ coeffs = self.coeffs m = hashlib.sha1() m.update(self.origin) m.update(self.order) for d in coeffs: m.update(str(d).encode("utf-8")) m.update(str(self.dx).encode("utf-8")) return m.hexdigest()[:keep_only]
[docs] def is_symbolic(self): """ Check if either the stencil leading factor or the stencil coefficients contains symbolic values. Returns ------- is_symbolic : bool """ return len(self.variables()) != 0
[docs] def is_centered(self, axe=None): """ Check if the stencil is centered (ie. stencil generate a centered finite differences scheme). Returns ------- is_centered : bool """ L, R = self.L, self.R if axe == None: return (L == R).all() else: return (L[axe] == R[axe]).all()
[docs] def is_cross(self): """ Check if the stencil is a cross (zeros everywhere except on a n-dimensional axis aligned cross). Examples: oo+oo o+oo oo+oo ++++ +++++ o+oo oo+oo o+oo Returns ------- is_cross : bool """ mask = np.ones(self.shape, dtype=bool) if self.dim == 1: return True else: access = self.origin.tolist() for i in range(self.dim): acc = [x for x in access] acc[i] = slice(0, self.shape[i]) mask[acc] = False return not np.any((self.coeffs != 0) * mask)
def __str__(self): return """ == {}D Stencil == origin: {} order: {} accuracy: {} error: {} shape: {} centered: {} cross: {} symbolic: {} variables: {} factor: {} coeffs: {} ================= """.format( self.dim, self.origin, self.order, self.accuracy(), self.error, self.shape, self.is_centered(), self.is_cross(), self.is_symbolic(), self.variables(), self.factor, self.coeffs, )
[docs] class CenteredStencil(Stencil): """ Centered stencil used for finite abitrary dimension differences. """ def __init__( self, coeffs, origin, order, dx=sm.Symbol("dx"), factor=1, error=None, **kwds ): """ Centered stencil used for finite abitrary dimension differences. See also -------- :class:`CenteredStencilGenerator`: Centered stencil generator. """ shape = np.asarray(coeffs.shape) if (shape % 2 == 0).any(): raise ValueError("Shape compnonent even!") if (origin != (shape - 1) // 2).any(): print(origin) print((shape - 1) // 2) raise ValueError("Origin is not centered!") super().__init__(coeffs, origin, order, dx, factor, error, **kwds)
[docs] def is_centered(self): return True
[docs] @staticmethod def from_stencil(stencil): s = stencil if not s.is_centered(): raise ValueError("Given stencil is not centered!") return CenteredStencil(s.coeffs, s.origin, s.order, s.dx, s.factor, s.error)